C     ALGORITHM 600, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2,
C     JUN., 1983, P. 258-259.
      SUBROUTINE QUINAT(N, X, Y, B, C, D, E, F)
C
      INTEGER N
      REAL X(N), Y(N), B(N), C(N), D(N), E(N), F(N)
C
C
C
C     QUINAT COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL QUINTIC SPLI
C     S(X) WITH KNOTS X(I) INTERPOLATING THERE TO GIVEN FUNCTION VALUES:
C               S(X(I)) = Y(I)  FOR I = 1,2, ..., N.
C     IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS A
C     POLYNOMIAL OF FIFTH DEGREE:
C     S(XX) = ((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I)    (*)
C           = ((((-F(I)*Q+E(I+1))*Q-D(I+1))*Q+C(I+1))*Q-B(I+1))*Q+Y(I+1)
C     WHERE  P = XX - X(I)  AND  Q = X(I+1) - XX.
C     (NOTE THE FIRST SUBSCRIPT IN THE SECOND EXPRESSION.)
C     THE DIFFERENT POLYNOMIALS ARE PIECED TOGETHER SO THAT S(X) AND
C     ITS DERIVATIVES UP TO S"" ARE CONTINUOUS.
C
C        INPUT:
C
C     N          NUMBER OF DATA POINTS, (AT LEAST THREE, I.E. N > 2)
C     X(1:N)     THE STRICTLY INCREASING OR DECREASING SEQUENCE OF
C                KNOTS.  THE SPACING MUST BE SUCH THAT THE FIFTH POWER
C                OF X(I+1) - X(I) CAN BE FORMED WITHOUT OVERFLOW OR
C                UNDERFLOW OF EXPONENTS.
C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS.
C
C        OUTPUT:
C
C     B,C,D,E,F  THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
C         (1:N)  SPECIFICALLY
C                B(I) = S'(X(I)), C(I) = S"(X(I))/2, D(I) = S"'(X(I))/6,
C                E(I) = S""(X(I))/24,  F(I) = S""'(X(I))/120.
C                F(N) IS NEITHER USED NOR ALTERED.  THE FIVE ARRAYS
C                B,C,D,E,F MUST ALWAYS BE DISTINCT.
C
C        OPTION:
C
C     IT IS POSSIBLE TO SPECIFY VALUES FOR THE FIRST AND SECOND
C     DERIVATIVES OF THE SPLINE FUNCTION AT ARBITRARILY MANY KNOTS.
C     THIS IS DONE BY RELAXING THE REQUIREMENT THAT THE SEQUENCE OF
C     KNOTS BE STRICTLY INCREASING OR DECREASING.  SPECIFICALLY:
C
C     IF X(J) = X(J+1) THEN S(X(J)) = Y(J) AND S'(X(J)) = Y(J+1),
C     IF X(J) = X(J+1) = X(J+2) THEN IN ADDITION S"(X(J)) = Y(J+2).
C
C     NOTE THAT S""(X) IS DISCONTINUOUS AT A DOUBLE KNOT AND, IN
C     ADDITION, S"'(X) IS DISCONTINUOUS AT A TRIPLE KNOT.  THE
C     SUBROUTINE ASSIGNS Y(I) TO Y(I+1) IN THESE CASES AND ALSO TO
C     Y(I+2) AT A TRIPLE KNOT.  THE REPRESENTATION (*) REMAINS
C     VALID IN EACH OPEN INTERVAL (X(I),X(I+1)).  AT A DOUBLE KNOT,
C     X(J) = X(J+1), THE OUTPUT COEFFICIENTS HAVE THE FOLLOWING VALUES:
C       Y(J) = S(X(J))          = Y(J+1)
C       B(J) = S'(X(J))         = B(J+1)
C       C(J) = S"(X(J))/2       = C(J+1)
C       D(J) = S"'(X(J))/6      = D(J+1)
C       E(J) = S""(X(J)-0)/24     E(J+1) = S""(X(J)+0)/24
C       F(J) = S""'(X(J)-0)/120   F(J+1) = S""'(X(J)+0)/120
C     AT A TRIPLE KNOT, X(J) = X(J+1) = X(J+2), THE OUTPUT
C     COEFFICIENTS HAVE THE FOLLOWING VALUES:
C       Y(J) = S(X(J))         = Y(J+1)    = Y(J+2)
C       B(J) = S'(X(J))        = B(J+1)    = B(J+2)
C       C(J) = S"(X(J))/2      = C(J+1)    = C(J+2)
C       D(J) = S"'((X(J)-0)/6    D(J+1) = 0  D(J+2) = S"'(X(J)+0)/6
C       E(J) = S""(X(J)-0)/24    E(J+1) = 0  E(J+2) = S""(X(J)+0)/24
C       F(J) = S""'(X(J)-0)/120  F(J+1) = 0  F(J+2) = S""'(X(J)+0)/120
C
      INTEGER I, M
      REAL B1, P, PQ, PQQR, PR, P2, P3, Q, QR, Q2, Q3, R, R2, S, T, U, V
C
      IF (N.LE.2) GO TO 190
C
C     COEFFICIENTS OF A POSITIVE DEFINITE, PENTADIAGONAL MATRIX,
C     STORED IN D,E,F FROM 2 TO N-2.
C
      M = N - 2
      Q = X(2) - X(1)
      R = X(3) - X(2)
      Q2 = Q*Q
      R2 = R*R
      QR = Q + R
      D(1) = 0.
      E(1) = 0.
      D(2) = 0.
      IF (Q.NE.0.) D(2) = 6.*Q*Q2/(QR*QR)
C
      IF (M.LT.2) GO TO 40
      DO 30 I=2,M
        P = Q
        Q = R
        R = X(I+2) - X(I+1)
        P2 = Q2
        Q2 = R2
        R2 = R*R
        PQ = QR
        QR = Q + R
        IF (Q) 20, 10, 20
   10   D(I+1) = 0.
        E(I) = 0.
        F(I-1) = 0.
        GO TO 30
   20   Q3 = Q2*Q
        PR = P*R
        PQQR = PQ*QR
        D(I+1) = 6.*Q3/(QR*QR)
        D(I) = D(I) + (Q+Q)*(15.*PR*PR+(P+R)*Q*(20.*PR+7.*Q2)+Q2*(8.*
     *   (P2+R2)+21.*PR+Q2+Q2))/(PQQR*PQQR)
        D(I-1) = D(I-1) + 6.*Q3/(PQ*PQ)
        E(I) = Q2*(P*QR+3.*PQ*(QR+R+R))/(PQQR*QR)
        E(I-1) = E(I-1) + Q2*(R*PQ+3.*QR*(PQ+P+P))/(PQQR*PQ)
        F(I-1) = Q3/PQQR
   30 CONTINUE
C
   40 IF (R.NE.0.) D(M) = D(M) + 6.*R*R2/(QR*QR)
C
C     FIRST AND SECOND ORDER DIVIDED DIFFERENCES OF THE GIVEN FUNCTION
C     VALUES, STORED IN B FROM 2 TO N AND IN C FROM 3 TO N
C     RESPECTIVELY. CARE IS TAKEN OF DOUBLE AND TRIPLE KNOTS.
C
      DO 60 I=2,N
        IF (X(I).NE.X(I-1)) GO TO 50
        B(I) = Y(I)
        Y(I) = Y(I-1)
        GO TO 60
   50   B(I) = (Y(I)-Y(I-1))/(X(I)-X(I-1))
   60 CONTINUE
      DO 80 I=3,N
        IF (X(I).NE.X(I-2)) GO TO 70
        C(I) = B(I)*0.5
        B(I) = B(I-1)
        GO TO 80
   70   C(I) = (B(I)-B(I-1))/(X(I)-X(I-2))
   80 CONTINUE
C
C     SOLVE THE LINEAR SYSTEM WITH C(I+2) - C(I+1) AS RIGHT-HAND SIDE.
C
      IF (M.LT.2) GO TO 100
      P = 0.
      C(1) = 0.
      E(M) = 0.
      F(1) = 0.
      F(M-1) = 0.
      F(M) = 0.
      C(2) = C(4) - C(3)
      D(2) = 1./D(2)
C
      IF (M.LT.3) GO TO 100
      DO 90 I=3,M
        Q = D(I-1)*E(I-1)
        D(I) = 1./(D(I)-P*F(I-2)-Q*E(I-1))
        E(I) = E(I) - Q*F(I-1)
        C(I) = C(I+2) - C(I+1) - P*C(I-2) - Q*C(I-1)
        P = D(I-1)*F(I-1)
   90 CONTINUE
C
  100 I = N - 1
      C(N-1) = 0.
      C(N) = 0.
      IF (N.LT.4) GO TO 120
      DO 110 M=4,N
C        I = N-2, ..., 2
        I = I - 1
        C(I) = (C(I)-E(I)*C(I+1)-F(I)*C(I+2))*D(I)
  110 CONTINUE
C
C     INTEGRATE THE THIRD DERIVATIVE OF S(X).
C
  120 M = N - 1
      Q = X(2) - X(1)
      R = X(3) - X(2)
      B1 = B(2)
      Q3 = Q*Q*Q
      QR = Q + R
      IF (QR) 140, 130, 140
  130 V = 0.
      T = 0.
      GO TO 150
  140 V = C(2)/QR
      T = V
  150 F(1) = 0.
      IF (Q.NE.0.) F(1) = V/Q
      DO 180 I=2,M
        P = Q
        Q = R
        R = 0.
        IF (I.NE.M) R = X(I+2) - X(I+1)
        P3 = Q3
        Q3 = Q*Q*Q
        PQ = QR
        QR = Q + R
        S = T
        T = 0.
        IF (QR.NE.0.) T = (C(I+1)-C(I))/QR
        U = V
        V = T - S
        IF (PQ) 170, 160, 170
  160   C(I) = C(I-1)
        D(I) = 0.
        E(I) = 0.
        F(I) = 0.
        GO TO 180
  170   F(I) = F(I-1)
        IF (Q.NE.0.) F(I) = V/Q
        E(I) = 5.*S
        D(I) = 10.*(C(I)-Q*S)
        C(I) = D(I)*(P-Q) + (B(I+1)-B(I)+(U-E(I))*P3-(V+E(I))*Q3)/PQ
        B(I) = (P*(B(I+1)-V*Q3)+Q*(B(I)-U*P3))/PQ -
     *   P*Q*(D(I)+E(I)*(Q-P))
  180 CONTINUE
C
C     END POINTS X(1) AND X(N).
C
      P = X(2) - X(1)
      S = F(1)*P*P*P
      E(1) = 0.
      D(1) = 0.
      C(1) = C(2) - 10.*S
      B(1) = B1 - (C(1)+S)*P
C
      Q = X(N) - X(N-1)
      T = F(N-1)*Q*Q*Q
      E(N) = 0.
      D(N) = 0.
      C(N) = C(N-1) + 10.*T
      B(N) = B(N) + (C(N)-T)*Q
  190 RETURN
      END
C
C
      SUBROUTINE QUINEQ(N, Y, B, C, D, E, F)
C
      INTEGER N
      REAL Y(N), B(N), C(N), D(N), E(N), F(N)
C
C
C
C     QUINEQ COMPUTES THE COEFFICIENTS OF QUINTIC NATURAL QUINTIC SPLINE
C     S(X) WITH EQUIDISTANT KNOTS X(I) INTERPOLATING THERE TO GIVEN
C     FUNCTION VALUES:
C               S(X(I)) = Y(I)  FOR I = 1,2, ..., N.
C     IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS
C     A POLYNOMIAL OF FIFTH DEGREE:
C     S(XX)=((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I)    (*)
C          =((((-F(I)*Q+E(I+1))*Q-D(I+1))*Q+C(I+1))*Q-B(I+1))*Q+Y(I+1)
C     WHERE  P = (XX - X(I))/X(I+1) - X(I))
C     AND    Q = (X(I+1) - XX)/(X(I+1) - X(I)).
C     (NOTE THE FIRST SUBSCRIPT IN THE SECOND EXPRESSION.)
C     THE DIFFERENT POLYNOMIALS ARE PIECED TOGETHER SO THAT S(X) AND
C     ITS DERIVATIVES UP TO S"" ARE CONTINUOUS.
C
C        INPUT:
C
C     N          NUMBER OF DATA POINTS, (AT LEAST THREE, I.E. N > 2)
C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS
C
C        OUTPUT:
C
C     B,C,D,E,F  THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
C         (1:N)  IF X(I+1) - X(I) = 1., THEN SPECIFICALLY:
C                B(I) = S'X(I)), C(I) = S"(X(I))/2, D(I) = S"'(X(I))/6,
C                E(I) = S""(X(I))/24,  F(I) = S""'(X(I)+0)/120.
C                F(N) IS NEITHER USED NOR ALTERED.  THE ARRAYS
C                Y,B,C,D MUST ALWAYS BE DISTINCT.  IF E AND F ARE
C                NOT WANTED, THE CALL QUINEQ(N,Y,B,C,D,D,D) MAY
C                BE USED TO SAVE STORAGE LOCATIONS.
C
      INTEGER I, M
      REAL P, Q, R, S, T, U, V
C
      IF (N.LE.2) GO TO 50
C
      M = N - 3
      P = 0.
      Q = 0.
      R = 0.
      S = 0.
      T = 0.
      D(M+1) = 0.
      D(M+2) = 0.
      IF (M.LE.0) GO TO 30
      DO 10 I=1,M
        U = P*R
        B(I) = 1./(66.-U*R-Q)
        R = 26. - U
        C(I) = R
        D(I) = Y(I+3) - 3.*(Y(I+2)-Y(I+1)) - Y(I) - U*S - Q*T
        Q = P
        P = B(I)
        T = S
        S = D(I)
   10 CONTINUE
C
      I = N - 2
      DO 20 M=4,N
C        I    = N-3, ..., 1
        I = I - 1
        D(I) = (D(I)-C(I)*D(I+1)-D(I+2))*B(I)
   20 CONTINUE
C
   30 M = N - 1
      Q = 0.
      R = D(1)
      T = R
      V = R
      DO 40 I=2,M
        P = Q
        Q = R
        R = D(I)
        S = T
        T = P - Q - Q + R
        F(I) = T
        U = 5.*(-P+Q)
        E(I) = U
        D(I) = 10.*(P+Q)
        C(I) = 0.5*(Y(I+1)+Y(I-1)+S-T) - Y(I) - U
        B(I) = 0.5*(Y(I+1)-Y(I-1)-S-T) - D(I)
   40 CONTINUE
C
      F(1) = V
      E(1) = 0.
      D(1) = 0.
      C(1) = C(2) - 10.*V
      B(1) = Y(2) - Y(1) - C(1) - V
      E(N) = 0.
      D(N) = 0.
      C(N) = C(N-1) + 10.*T
      B(N) = Y(N) - Y(N-1) + C(N) - T
   50 RETURN
      END
C
C
C
C
      SUBROUTINE QUINDF(N, X, Y, B, C, D, E, F)
C
      INTEGER N
      REAL X(N), Y(N), B(N), C(N), D(N), E(N), F(N)
C
C
C
C     QUINDF COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL QUINTIC
C     SPLINE S(X) WITH KNOTS X(I) FOR WHICH THE FUNCTION VALUES Y(I) AND
C     THE FIRST DERIVATIVES B(I) ARE SPECIFIED AT X(I), I = 1,2, ...,N.
C     IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS
C     A POLYNOMIAL OF FIFTH DEGREE:
C     S(XX)=((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I)    (*)
C     WHERE  P = XX - X(I).
C
C        INPUT:
C
C     N          NUMBER OF DATA POINTS, (AT LEAST TWO, I.E. N > 1)
C     X(1:N)     THE STRICTLY INCREASING OR DECREASING SEQUENCE OF
C                KNOTS.  THE SPACING MUST BE SUCH THAT THE FIFTH
C                POWER OF X(I+1) - X(I) CAN BE FORMED WITHOUT
C                OVERFLOW OR UNDERFLOW OF EXPONENTS
C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS
C     B(1:N)     THE PRESCRIBED DERIVATIVE VALUES AT THE KNOTS
C
C        OUTPUT:
C
C     C,D,E,F    THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
C         (1:N)  E(N) AND F(N) ARG NEITHER USED NOR ALTERED.
C                THE ARRAYS C,D,E,F MUST ALWAYS BE DISTINCT.
C
      INTEGER I, M, N1
      REAL CC, G, H, HH, H2, P, PP, Q, QQ, R, RR
C
      IF (N.LE.1) GO TO 40
      N1 = N - 1
      CC = 0.
      HH = 0.
      PP = 0.
      QQ = 0.
      RR = 0.
      G = 0.
      DO 10 I=1,N1
        H = 1./(X(I+1)-X(I))
        H2 = H*H
        D(I) = 3.*(HH+H) - G*HH
        P = (Y(I+1)-Y(I))*H2*H
        Q = (B(I+1)+B(I))*H2
        R = (B(I+1)-B(I))*H2
        CC = 10.*(P-PP) - 5.*(Q-QQ) + R + RR + G*CC
        C(I) = CC
        G = H/D(I)
        HH = H
        PP = P
        QQ = Q
        RR = R
   10 CONTINUE
C
      C(N) = (-10.*PP+5.*QQ+RR+G*CC)/(3.*HH-G*HH)
      I = N
      DO 20 M=1,N1
C        I      = N-1, ..., 1
        I = I - 1
        D(I+1) = 1./(X(I+1)-X(I))
        C(I) = (C(I)+C(I+1)*D(I+1))/D(I)
   20 CONTINUE
C
      DO 30 I=1,N1
        H = D(I+1)
        P = (((Y(I+1)-Y(I))*H-B(I))*H-C(I))*H
        Q = ((B(I+1)-B(I))*H-C(I)-C(I))*H
        R = (C(I+1)-C(I))*H
        G = Q - 3.*P
        RR = R - 3.*(P+G)
        QQ = -RR - RR + G
        F(I) = RR*H*H
        E(I) = QQ*H
        D(I) = -RR - QQ + P
   30 CONTINUE
C
      D(N) = 0.
      E(N) = 0.
      F(N) = 0.
   40 RETURN
      END
C    DRIVER PROGRAM FOR TEST OF QUINAT, QUINEQ AND QUINDF
C       FOLLOWS.
C
      INTEGER N,NM1,M,MM,MM1,I,K,J,JJ
      REAL Z
      REAL X(200),Y(200),B(200),BB(200),CC(200),DD(200),EE(200),
     *                 FF(200),A(200,6),C(6),DIFF(5),COM(5)
C
C     N          NUMBER OF DATA POINTS.
C     M          2*M-1 IS ORDER OF SPLINE.
C                   M = 3 ALWAYS FOR QUINTIC SPLINE.
C     NN,NM1,MM,
C     MM1,I,K,
C     J,JJ       TEMPORARY INTEGER VARIABLES.
C     Z,P        TEMPORARY REAL VARIABLES.
C     X(1:N)     THE SEQUENCE OF KNOTS.
C     Y(1:N)     THE PRESCRIBED FUNCTION VALUES AT THE KNOTS.
C     B(1:N)     THE PRESCRIBED DERIVATIVE VALUES AT THE KNOTS.
C     BB,CC,DD,
C     EE,FF(1:N) THE COMPUTED SPLINE COEFFICIENTS
C     A(1:N,1:6) TWO DIMENSIONAL ARRAY WHOSE COLUMNS ARE
C                   Y, BB, CC, DD, EE, FF.
C     DIFF(1:5)  MAXIMUM VALUES OF DIFFERENCES OF VALUES AND
C                   DERIVATIVES TO RIGHT AND LEFT OF KNOTS.
C     COM(1:5)   MAXIMUM VALUES OF COEFFICIENTS.
C
C
C     TEST OF QUINAT WITH NONEQUIDISTANT KNOTS AND
C        EQUIDISTANT KNOTS FOLLOWS.
C
      NOUT=6
      WRITE(NOUT,10)
   10 FORMAT(50H1         TEST OF QUINAT WITH NONEQUIDISTANT KNOTS)
      N = 5
      X(1) = -3.0
      X(2) = -1.0
      X(3) =  0.0
      X(4) =  3.0
      X(5) =  4.0
      Y(1) =  7.0
      Y(2) = 11.0
      Y(3) = 26.0
      Y(4) = 56.0
      Y(5) = 29.0
      M = 3
      MM = 2*M
      MM1 = MM - 1
      WRITE(NOUT,15) N,M
   15 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
      DO 20 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
   20 CONTINUE
      DO 30 I = 1,MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
   30 CONTINUE
      DO 70 K = 1,N
         DO 35 I = 1,MM
            C(I) = A(K,I)
   35 CONTINUE
      WRITE(NOUT,40) K
   40 FORMAT(40H ---------------------------------------,I3,
     *  45H --------------------------------------------)
      WRITE(NOUT,45) X(K)
   45 FORMAT(F12.8)
      IF (K .EQ. N) WRITE(NOUT,55) C(1)
      IF (K .EQ. N) GO TO 75
      WRITE(NOUT,55) (C(I),I=1,MM)
      DO 50 I = 1,MM1
         IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
   50 CONTINUE
   55 FORMAT(6F16.8)
      Z = X(K+1) - X(K)
      DO 60 I = 2,MM
         DO 60 JJ = I,MM
            J = MM + I - JJ
            C(J-1) = C(J)*Z + C(J-1)
   60 CONTINUE
      WRITE(NOUT,55) (C(I),I=1,MM)
      DO 65 I = 1,MM1
         IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 65
         Z = ABS(C(I) - A(K+1,I))
         IF (Z .GT. DIFF(I)) DIFF(I) = Z
   65 CONTINUE
   70 CONTINUE
   75 WRITE(NOUT,80)
   80 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,85) (DIFF(I),I=1,MM1)
   85 FORMAT(5E18.9)
      WRITE(NOUT,90)
   90 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
      WRITE(NOUT,95) (COM(I),I=1,MM1)
   95 FORMAT(5F16.8)
      M = 3
      DO 200 N = 10,100,10
      MM = 2*M
      MM1 = MM - 1
      NM1 = N -1
      DO 100 I = 1,NM1,2
         X(I)   = I
         X(I+1) = I + 1
         Y(I)   = 1.
         Y(I+1) = 0.
  100 CONTINUE
      IF (MOD(N,2) .EQ. 0) GOTO 105
      X(N) = N
      Y(N) = 1.
  105 CONTINUE
      WRITE(NOUT,110) N,M
  110 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
      DO 115 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  115 CONTINUE
      DO 120 I = 1, MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  120 CONTINUE
      DO 165 K = 1,N
         DO 125 I = 1,MM
            C(I) = A(K,I)
  125    CONTINUE
         IF (N .GT. 10) GOTO 140
         WRITE(NOUT,130) K
  130 FORMAT(40H ---------------------------------------,I3,
     *       45H --------------------------------------------)
         WRITE(NOUT,135) X(K)
  135 FORMAT(F12.8)
         IF (K .EQ. N) WRITE(NOUT,150) C(1)
  140    CONTINUE
         IF (K .EQ. N) GO TO 170
         IF (N .LE. 10) WRITE(NOUT,150) (C(I), I=1,MM)
         DO 145 I = 1,MM1
            IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  145    CONTINUE
  150 FORMAT(6F16.8)
         Z = X(K+1) - X(K)
         DO 155 I = 2,MM
            DO 155 JJ = I,MM
               J = MM + I - JJ
               C(J-1) = C(J)*Z + C(J-1)
  155    CONTINUE
         IF (N .LE. 10) WRITE(NOUT,150) (C(I), I=1,MM)
         DO 160 I = 1,MM1
            IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 160
            Z = ABS(C(I) - A(K+1,I))
            IF (Z .GT. DIFF(I)) DIFF(I) = Z
  160    CONTINUE
  165 CONTINUE
  170 WRITE(NOUT,175)
  175 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,180) (DIFF(I),I=1,MM1)
  180 FORMAT(5E18.9)
      WRITE(NOUT,185)
  185 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
      WRITE(NOUT,190) (COM(I),I=1,MM1)
  190 FORMAT(5E18.9)
  200 CONTINUE
C
C
C     TEST OF QUINEQ FOLLOWS.
C
      WRITE(NOUT,210)
  210 FORMAT(18H1   TEST OF QUINEQ)
      M = 3
      DO 400 N = 10,100,10
      MM = 2*M
      MM1 = MM - 1
      NM1 = N -1
      DO 300 I = 1,NM1,2
         X(I)   = I
         X(I+1) = I + 1
         Y(I)   = 1.
         Y(I+1) = 0.
  300 CONTINUE
      IF (MOD(N,2) .EQ. 0) GOTO 305
      X(N) = FLOAT(N)
      Y(N) = 1.
  305 CONTINUE
      WRITE(NOUT,310) N,M
  310 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINEQ(N,Y,BB,CC,DD,EE,FF)
      DO 315 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  315 CONTINUE
      DO 320 I = 1, MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  320 CONTINUE
      DO 365 K = 1,N
         DO 325 I = 1,MM
            C(I) = A(K,I)
  325    CONTINUE
         IF (N .GT. 10) GOTO 340
         WRITE(NOUT,330) K
  330 FORMAT(40H ---------------------------------------,I3,
     *       45H --------------------------------------------)
         WRITE(NOUT,335) X(K)
  335 FORMAT(F12.8)
         IF (K .EQ. N) WRITE(NOUT,350) C(1)
  340    CONTINUE
         IF (K .EQ. N) GO TO 370
         IF (N .LE. 10) WRITE(NOUT,350) (C(I), I=1,MM)
         DO 345 I = 1,MM1
            IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  345    CONTINUE
  350 FORMAT(6F16.8)
         Z = 1.
         DO 355 I = 2,MM
            DO 355 JJ = I,MM
               J = MM + I - JJ
               C(J-1) = C(J)*Z + C(J-1)
  355    CONTINUE
         IF (N .LE. 10) WRITE(NOUT,350) (C(I), I=1,MM)
         DO 360 I = 1,MM1
            IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 360
            Z = ABS(C(I) - A(K+1,I))
            IF (Z .GT. DIFF(I)) DIFF(I) = Z
  360    CONTINUE
  365 CONTINUE
  370 WRITE(NOUT,375)
  375 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,380) (DIFF(I),I=1,MM1)
  380 FORMAT(5E18.9)
      WRITE(NOUT,385)
  385 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
      WRITE(NOUT,390) (COM(I),I=1,MM1)
  390 FORMAT(5E18.9)
  400 CONTINUE
C
C
C     TEST OF QUINDF WITH NONEQUIDISTANT KNOTS FOLLOWS.
C
      WRITE(NOUT,410)
  410 FORMAT(50H1         TEST OF QUINDF WITH NONEQUIDISTANT KNOTS)
      N = 5
      X(1) = -3.0
      X(2) = -1.0
      X(3) =  0.0
      X(4) =  3.0
      X(5) =  4.0
      Y(1) =  7.0
      Y(2) = 11.0
      Y(3) = 26.0
      Y(4) = 56.0
      Y(5) = 29.0
      B(1) =  2.0
      B(2) = 15.0
      B(3) = 10.0
      B(4) =-27.0
      B(5) =-30.0
      M = 3
      MM = 2*M
      MM1 = MM - 1
      WRITE(NOUT,415) N,M
  415 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINDF(N,X,Y,B,CC,DD,EE,FF)
      DO 420 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = B(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  420 CONTINUE
      DO 430 I = 1,MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  430 CONTINUE
      DO 470 K = 1,N
         DO 435 I = 1,MM
            C(I) = A(K,I)
  435 CONTINUE
      WRITE(NOUT,440) K
  440 FORMAT(40H ---------------------------------------,I3,
     *  45H --------------------------------------------)
      WRITE(NOUT,445) X(K)
  445 FORMAT(F12.8)
      IF (K .EQ. N) WRITE(NOUT,455) C(1)
      IF (K .EQ. N) GO TO 475
      WRITE(NOUT,455) (C(I),I=1,MM)
      DO 450 I = 1,MM1
         IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  450 CONTINUE
  455 FORMAT(6F16.8)
      Z = X(K+1) - X(K)
      DO 460 I = 2,MM
         DO 460 JJ = I,MM
            J = MM + I - JJ
            C(J-1) = C(J)*Z + C(J-1)
  460 CONTINUE
      WRITE(NOUT,455) (C(I),I=1,MM)
      DO 465 I = 1,MM1
         IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 465
         Z = ABS(C(I) - A(K+1,I))
         IF (Z .GT. DIFF(I)) DIFF(I) = Z
  465 CONTINUE
  470 CONTINUE
  475 WRITE(NOUT,480)
  480 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,485) (DIFF(I),I=1,MM1)
  485 FORMAT(5E18.9)
      WRITE(NOUT,490)
  490 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
      WRITE(NOUT,495) (COM(I),I=1,MM1)
  495 FORMAT(5F16.8)
      M = 3
      DO 600 N = 10,100,10
      MM = 2*M
      MM1 = MM - 1
      P = 0.0
      DO 500 I = 1,N
         P = P + ABS(SIN((FLOAT(I)))) + 0.001*FLOAT(I)
         X(I) = P
         Y(I) = COS((FLOAT(I))) - 0.5
         B(I) = COS((FLOAT(2*I))) - 0.5
  500 CONTINUE
      WRITE(NOUT,510) N,M
  510 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINDF(N,X,Y,B,CC,DD,EE,FF)
      DO 515 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = B(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  515 CONTINUE
      DO 520 I = 1, MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  520 CONTINUE
      DO 565 K = 1,N
         DO 525 I = 1,MM
            C(I) = A(K,I)
  525    CONTINUE
         IF (N .GT. 10) GOTO 540
         WRITE(NOUT,530) K
  530 FORMAT(40H ---------------------------------------,I3,
     *       45H --------------------------------------------)
         WRITE(NOUT,535) X(K)
  535 FORMAT(F12.8)
         IF (K .EQ. N) WRITE(NOUT,550) C(1)
  540    CONTINUE
         IF (K .EQ. N) GO TO 570
         IF (N .LE. 10) WRITE(NOUT,550) (C(I), I=1,MM)
         DO 545 I = 1,MM1
            IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  545    CONTINUE
  550 FORMAT(6F16.8)
         Z = X(K+1) - X(K)
         DO 555 I = 2,MM
            DO 555 JJ = I,MM
               J = MM + I - JJ
               C(J-1) = C(J)*Z + C(J-1)
  555    CONTINUE
         IF (N .LE. 10) WRITE(NOUT,550) (C(I), I=1,MM)
         DO 560 I = 1,MM1
            IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 560
            Z = ABS(C(I) - A(K+1,I))
            IF (Z .GT. DIFF(I)) DIFF(I) = Z
  560    CONTINUE
  565 CONTINUE
  570 WRITE(NOUT,575)
  575 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,580) (DIFF(I),I=1,MM1)
  580 FORMAT(5E18.9)
      WRITE(NOUT,585)
  585 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
      WRITE(NOUT,590) (COM(I),I=1,MM1)
  590 FORMAT(5E18.9)
  600 CONTINUE
C
C
C     TEST OF QUINAT WITH NONEQUIDISTANT DOUBLE KNOTS FOLLOWS.
C
      WRITE(NOUT,610)
  610 FORMAT(50H1  TEST OF QUINAT WITH NONEQUIDISTANT DOUBLE KNOTS)
      N = 5
      X(1) = -3.
      X(2) = -3.
      X(3) = -1.
      X(4) = -1.
      X(5) = 0.
      X(6) = 0.
      X(7) = 3.
      X(8) = 3.
      X(9) = 4.
      X(10) = 4.
      Y(1) = 7.
      Y(2) = 2.
      Y(3) = 11.
      Y(4) = 15.
      Y(5) = 26.
      Y(6) = 10.
      Y(7) = 56.
      Y(8) = -27.
      Y(9) = 29.
      Y(10) = -30.
      M = 3
      NN = 2*N
      MM = 2*M
      MM1 = MM - 1
      WRITE(NOUT,615) N,M
  615 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINAT(NN,X,Y,BB,CC,DD,EE,FF)
      DO 620 I = 1,NN
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  620 CONTINUE
      DO 630 I = 1,MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  630 CONTINUE
      DO 670 K = 1,NN
         DO 635 I = 1,MM
            C(I) = A(K,I)
  635 CONTINUE
      WRITE(NOUT,640) K
  640 FORMAT(40H ---------------------------------------,I3,
     *  45H --------------------------------------------)
      WRITE(NOUT,645) X(K)
  645 FORMAT(F12.8)
      IF (K .EQ. NN) WRITE(NOUT,655) C(1)
      IF (K .EQ. NN) GO TO 675
      WRITE(NOUT,655) (C(I),I=1,MM)
      DO 650 I = 1,MM1
         IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  650 CONTINUE
  655 FORMAT(6F16.8)
      Z = X(K+1) - X(K)
      DO 660 I = 2,MM
         DO 660 JJ = I,MM
            J = MM + I - JJ
            C(J-1) = C(J)*Z + C(J-1)
  660 CONTINUE
      WRITE(NOUT,655) (C(I),I=1,MM)
      DO 665 I = 1,MM1
         IF ((K .GE. NN-1) .AND. (I .NE. 1)) GO TO 665
         Z = ABS(C(I) - A(K+1,I))
         IF (Z .GT. DIFF(I)) DIFF(I) = Z
  665 CONTINUE
  670 CONTINUE
  675 WRITE(NOUT,680)
  680 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,685) (DIFF(I),I=1,MM1)
  685 FORMAT(5E18.9)
      WRITE(NOUT,690)
  690 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
      WRITE(NOUT,695) (COM(I),I=1,MM1)
  695 FORMAT(5F16.8)
      M = 3
      DO 800 N = 10,100,10
      NN = 2*N
      MM = 2*M
      MM1 = MM - 1
      P = 0.0
      DO 700 I = 1,N
         P = P + ABS(SIN(FLOAT(I)))
         X(2*I-1) = P
         X(2*I)   = P
         Y(2*I-1) = COS(FLOAT(I)) - 0.5
         Y(2*I)   = COS(FLOAT(2*I)) - 0.5
  700 CONTINUE
      WRITE(NOUT,710) N,M
  710 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINAT(NN,X,Y,BB,CC,DD,EE,FF)
      DO 715 I = 1,NN
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  715 CONTINUE
      DO 720 I = 1, MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  720 CONTINUE
      DO 765 K = 1,NN
         DO 725 I = 1,MM
            C(I) = A(K,I)
  725    CONTINUE
         IF (N .GT. 10)  GOTO 740
         WRITE(NOUT,730) K
  730 FORMAT(40H ---------------------------------------,I3,
     *       45H --------------------------------------------)
         WRITE(NOUT,735) X(K)
  735 FORMAT(F12.8)
         IF (K .EQ. NN) WRITE(NOUT,750) C(1)
  740    CONTINUE
         IF (K .EQ. NN) GO TO 770
         IF (N .LE. 10) WRITE(NOUT,750) (C(I), I=1,MM)
         DO 745 I = 1,MM1
            IF (ABS(A(K,I)) .GT.  COM(I)) COM(I) = ABS(A(K,I))
  745    CONTINUE
  750 FORMAT(6F16.8)
         Z = X(K+1) - X(K)
         DO 755 I = 2,MM
            DO 755 JJ = I,MM
               J = MM + I - JJ
               C(J-1) = C(J)*Z + C(J-1)
  755    CONTINUE
         IF (N .LE. 10) WRITE(NOUT,750) (C(I), I=1,MM)
         DO 760 I = 1,MM1
            IF ((K .GE. NN-1) .AND. (I .NE. 1)) GO TO 760
            Z = ABS(C(I) - A(K+1,I))
            IF (Z .GT. DIFF(I)) DIFF(I) = Z
  760    CONTINUE
  765 CONTINUE
  770 WRITE(NOUT,775)
  775 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,780) (DIFF(I),I=1,MM1)
  780 FORMAT(5E18.9)
      WRITE(NOUT,785)
  785 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
      WRITE(NOUT,790) (COM(I),I=1,MM1)
  790 FORMAT(5E18.9)
  800 CONTINUE
C
C
C     TEST OF QUINAT WITH NONEQUIDISTANT KNOTS, ONE DOUBLE KNOT,
C        ONE TRIPLE KNOT, FOLLOWS.
C
      WRITE(NOUT,805)
      WRITE(NOUT,810)
  805 FORMAT(51H1         TEST OF QUINAT WITH NONEQUIDISTANT KNOTS,)
  810 FORMAT(40H             ONE DOUBLE, ONE TRIPLE KNOT)
      N = 8
      X(1) = -3.
      X(2) = -1.
      X(3) =  -1.
      X(4) =  0.
      X(5) =  3.
      X(6) = 3.
      X(7) = 3.
      X(8) = 4.
      Y(1) =  7.
      Y(2) = 11.
      Y(3) = 15.
      Y(4) = 26.
      Y(5) = 56.
      Y(6) = -30.
      Y(7) =  -7.
      Y(8) =  29.
      M = 3
      MM = 2*M
      MM1 = MM - 1
      WRITE(NOUT,815) N,M
  815 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
      DO 820 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  820 CONTINUE
      DO 830 I = 1,MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  830 CONTINUE
      DO 870 K = 1,N
         DO 835 I = 1,MM
            C(I) = A(K,I)
  835 CONTINUE
      WRITE(NOUT,840) K
  840 FORMAT(40H ---------------------------------------,I3,
     *  45H --------------------------------------------)
      WRITE(NOUT,845) X(K)
  845 FORMAT(F12.8)
      IF (K .EQ. N) WRITE(NOUT,855) C(1)
      IF (K .EQ. N) GO TO 875
      WRITE(NOUT,855) (C(I),I=1,MM)
      DO 850 I = 1,MM1
         IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  850 CONTINUE
  855 FORMAT(6F16.8)
      Z = X(K+1) - X(K)
      DO 860 I = 2,MM
         DO 860 JJ = I,MM
            J = MM + I - JJ
            C(J-1) = C(J)*Z + C(J-1)
  860 CONTINUE
      WRITE(NOUT,855) (C(I),I=1,MM)
      DO 865 I = 1,MM1
         IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 865
         Z = ABS(C(I) - A(K+1,I))
         IF (Z .GT. DIFF(I)) DIFF(I) = Z
  865 CONTINUE
  870 CONTINUE
  875 WRITE(NOUT,880)
  880 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,885) (DIFF(I),I=1,MM1)
  885 FORMAT(5E18.9)
      WRITE(NOUT,890)
  890 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
      WRITE(NOUT,895) (COM(I),I=1,MM1)
  895 FORMAT(5F16.8)
C
C
C     TEST OF QUINAT WITH NONEQUIDISTANT KNOTS, TWO DOUBLE KNOTS,
C        ONE TRIPLE KNOT,FOLLOWS.
C
      WRITE(NOUT,905)
      WRITE(NOUT,910)
  905 FORMAT(51H1         TEST OF QUINAT WITH NONEQUIDISTANT KNOTS,)
  910 FORMAT(40H             TWO DOUBLE, ONE TRIPLE KNOT)
      N = 10
      X(1) = 0.
      X(2) = 2.
      X(3) = 2.
      X(4) = 3.
      X(5) = 3.
      X(6) = 3.
      X(7) = 5.
      X(8) = 8.
      X(9) = 9.
      X(10)= 9.
      Y(1) = 163.
      Y(2) = 237.
      Y(3) = -127.
      Y(4) = 119.
      Y(5) = -65.
      Y(6) = 192.
      Y(7) = 293.
      Y(8) =  326.
      Y(9) = 0.
      Y(10)= -414.0
      M = 3
      MM = 2*M
      MM1 = MM - 1
      WRITE(NOUT,915) N,M
  915 FORMAT(5H-N = ,I3,7H    M =,I2)
      CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
      DO 920 I = 1,N
         A(I,1) = Y(I)
         A(I,2) = BB(I)
         A(I,3) = CC(I)
         A(I,4) = DD(I)
         A(I,5) = EE(I)
         A(I,6) = FF(I)
  920 CONTINUE
      DO 930 I = 1,MM1
         DIFF(I) = 0.0
         COM(I) = 0.0
  930 CONTINUE
      DO 970 K = 1,N
         DO 935 I = 1,MM
            C(I) = A(K,I)
  935 CONTINUE
      WRITE(NOUT,940) K
  940 FORMAT(40H ---------------------------------------,I3,
     *  45H --------------------------------------------)
      WRITE(NOUT,945) X(K)
  945 FORMAT(F12.8)
      IF (K .EQ. N) WRITE(NOUT,955) C(1)
      IF (K .EQ. N) GO TO 975
      WRITE(NOUT,955) (C(I),I=1,MM)
      DO 950 I = 1,MM1
         IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
  950 CONTINUE
  955 FORMAT(6F16.8)
      Z = X(K+1) - X(K)
      DO 960 I = 2,MM
         DO 960 JJ = I,MM
            J = MM + I - JJ
            C(J-1) = C(J)*Z + C(J-1)
  960 CONTINUE
      WRITE(NOUT,955) (C(I),I=1,MM)
      DO 965 I = 1,MM1
         IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 965
         Z = ABS(C(I) - A(K+1,I))
         IF (Z .GT. DIFF(I)) DIFF(I) = Z
  965 CONTINUE
  970 CONTINUE
  975 WRITE(NOUT,980)
  980 FORMAT(41H  MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
      WRITE(NOUT,985) (DIFF(I),I=1,MM1)
  985 FORMAT(5E18.9)
      WRITE(NOUT,990)
  990 FORMAT(42H  MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
      IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
      WRITE(NOUT,995) (COM(I),I=1,MM1)
  995 FORMAT(5F16.8)
      STOP
      END
